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We explore the collective density oscillations of a collection of charged massive Dirac particles, 
in one, two and three dimensions and their one dimensional superlattice. We calculate the long 
wavelength limit of the dynamical polarization function analytically, and use the random phase 
approximation to obtain the plasmon dispersion. The density dependence of the long wavelength 
plasmon frequency in massive Dirac systems is found to be different compared to systems with 
parabolic, and gapless Dirac dispersion. We also calculate the long wavelength plasmon dispersion 
of a Id metamaterial made from Id and 2d massive Dirac plasma. Our analytical results will be useful 
for exploring the use of massive Dirac materials as electrostatically tunable plasmonic metamaterials 
and can be experimentally verified by infrared spectroscopy as in the case of graphene [L. Ju. et. 
al, Nat. Nanotechnol. 6, 630 (2011)]. 
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I. INTRODUCTION 

The collective density oscillations of electrons liquids, 
i.e., plasmons, offer a powerful tool for exploring electron- 
electron interaction effects in various systems 1,2 and have 
also motivated several potential applications in opti¬ 
cal metamaterials, nanophotonic lasers and amplifiers, 
biochemical sensing, and antennas transmitting and re¬ 
ceiving light signals at the nanoscale 3,4 . The collective 
modes of ordinary (Schrodinger) electrons with parabolic 
dispersion 1,2,5 including spin-orbit coupling 6 and spin 
polarization' have been extensively studied in metals and 
doped semiconductors. Since the discovery of graphene, 
there has been huge interest in plasmons of Dirac materi¬ 
als and particularly in graphene 8-14 as it offers a tunable 
plasmon spectrum via electrostatic control of its carrier 
concentration, and higher plasmon lifetimes due to high 
mobility. 

There have been several studies on plasmons in gap¬ 
less two dimensional (2d) and three dimensional (3d) 
Dirac systems in the context of graphene 15-19 , topolog¬ 
ical insulators 20-22 , Weyl semimetals 23 , and in gapped 
2d massive Dirac systems 24 in the context of buckled 
honeycomb structures such as silicene 25-27 . In addi¬ 
tion to this plasmons in periodic arrays of parabolic 
systems 28-31 , and massless Dirac plasma layers 32-34 have 
also been investigated. A nretanraterial made up of pe¬ 
riodic graphene micro-ribbon arrays was used in Ref. [8] 
to demonstrate tunable terahertz plasmon excitations in 
graphene. Comparatively massive Dirac system in var¬ 
ious dimensions, and its multilayers/superlattice have 
been relatively less explored, and consequently are the 
subject of this article. 

In this article we study the plasmon frequency, and its 
density dependence for a massive Dirac plasma (MDP) 
interacting via the long ranged Coulomb interaction, in 
one-, two- and three dimensions. Additionally we also 
calculate the plasmon dispersion for metamaterials made 


of Id nano-ribbons and 2d layers of MDP. The gapless 
Dirac systems were studied in Ref. [18], which serves as 
a check of all of our calculations in the limit of the van¬ 
ishing band gap. We find that while the long-wavelength 
plasmon frequency in MDP is essentially quantum me¬ 
chanical in nature with l/\/h appearing explicitly in the 
plasmon dispersion as in the case of gapless Dirac plasma 
(GDP), the scaling of the plasmon frequency with den¬ 
sity is different for MDP, GDP and parabolic systems. 
Note that in systems with non-relativistic parabolic dis¬ 
persion, the long wavelength plasmon frequency is ‘classi¬ 
cal’ and quantum corrections (arising from the self energy 
and vertex corrections in the polarization function) show 
up only in higher order terms. The aim of this work 
is to illustrate the key differences between the density 
dependences of plasmon dispersions in one-, two- and 
three-dimensional systems which arise due to the rela¬ 
tivistic (Dirac) or non-relativistic (Schrodinger) nature 
of the electrons and due to the presence of a finite gap in 
MDP. 

This article is organized as follows: In Sec. II, we in¬ 
troduce the random phase approximation (RPA) ‘recipe’ 
for calculating the plasmon frequency and explicitly cal¬ 
culate the long wavelength limit of the dynamical polar¬ 
ization function for MDP, GDP and parabolic dispersion 
systems. This allows us to obtain and discuss similarities 
and differences in the long wavelength plasmon frequen¬ 
cies in Sec. III. Next we consider the plasmons arising 
in a periodic array of MDP nano-ribbons and layers in 
Sec. IV and compare the results with GDP and parabolic 
dispersion systems. Finally, in Sec. V, we summarize our 
findings. 


II. POLARIZATION FUNCTION 

Within the RPA, the collective plasmon modes of an 
electron system emerge as poles of the density density 
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response functions (also called polarization function or 
the Lindhard function) and coincide with the zeros of the 
complex longitudinal “dielectric function” e(q,ui), i.e. 


e(q, w) = 1 - v q Il(q, w) = 0 , (1) 


where v q is the Fourier transform of the Coulomb interac¬ 
tion, and II is the total non-interacting polarizability of 
the system. The Fourier transform of the Coulomb inter¬ 
action v{r ) = e 2 /(«x), in the appropriate d-dimensional 
space is given by 
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where k is the background material dependent dielec¬ 
tric constant, and K$ denotes the zeroth order modi¬ 
fied Bessel function of the second kind. Note that in 
one dimension, the length scale a characterizes the lat¬ 
eral confinement size (say radius of the Id ribbon), and 
v q « —2e 2 In (qa)/n for qa <C 1, while v q = e 2 /(nq 2 a 2 ) 
for qa 1. 

The polarization function for the massive Dirac mate¬ 
rial is given by 
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where k' = k + q, A, A' = ±1 denotes the con¬ 
duction (particle) and valence (hole) bands, E k = 
hv F \^/k 2 + ( A/Hv F ) 2 \ with 2A being the energy gap, 
n F (x) is the Fermi function and 2F) Vi- y(k, k') = 1 + 
AA'[k • k' + A 2 ]/(E k E k >) is the overlap function, with 
x = x/hv F . The factor g s (= 2) is the spin degeneracy 
factor and g v is the valley (or pseudo spin) degeneracy 
factor (e.g. g v = 2 for graphene and other Dirac materi¬ 
als with honeycomb lattice structure). Given the general 
relation II(g, —w) = II(g,u;)*, and the fact that the po¬ 
larization function depends only on the absolute value 
of the Fermi energy e F , we only present the results for 
> 0 and to > 0. Furthermore, we work at zero tem¬ 
perature so that the Fermi functions can be replaced by 
Heaviside step functions, i.e., n F (x) = 0(ef — x). 


Depending upon the placement of the Fermi energy 
s F . we can split our polarization function into two parts, 
namely the intrinsic (ep < A) and extrinsic (ep > A) 
polarization: 


where 

J d d k 0(D 2 - A 2 - k 2 ) (5) 
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Here the upper and lower signs correspond to intraband 
and interband electron-hole transitions respectively and 
the parameter D defines the integration limits via the 0 
function. Since we are interested in the long wavelength 
(q —> 0) plasmon dispersion, we evaluate Eq. (3) in the 
dynamical limit (q —> 0 first and then to —» 0), to lowest 
order in q 2 /cj 2 , just above the intra-band particle hole 
continuum. 

We mention at the outset that we will use the super¬ 
script (p), (g) and (m) to refer to systems with parabolic, 
gapless (or massless) Dirac, and massive Dirac systems, 
respectively. Note that the electronic density for any d- 
dimensional system, in terms of its Fermi wavevector is 
given by 
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where T(a;) is the Gamma function. However the Fermi 
wavevectors for parabolic, massive Dirac, and gapless 
Dirac systems are expressed differently in terms of the 
Fermi energy and are given by k F = y/2me F /h 1 k F = 
\/£ F — A 2 /Hv F and k F = e F /Hv F , respectively. 

For systems with parabolic dispersion (E k = 
h 2 k 2 /2m p ), Eq. (3) can be evaluated in the dynamical 
long wavelength limit, upto leading order in q just above 
the particle-hole continuum to obtain 
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For massive Dirac systems, E k = hv f \Jk 2 + A 2 (where 
A = A/ hv F ) in all dimensions, and the dynamical long 
wavelength limit of the Lindhard function, just above the 
particle-hole continuum, is given by 

n (m) _ 9s9vV f 7r d/2 _ ^f 9 2 n 

~ h(2n) d T(l + d/2)^^^ + 

( 8 ) 



For massless Dirac systems such as graphene, A —> 0, 
and E k = hv F k in all dimensions, and Eq. (8) reduces to 

gsgyvpk^r 1 n d / 2 (<f_\ 

H(2n) d T{l + d/2) w 2 \w 4 / 

(9) 



n (q,u>) = -x 00 (q , ,w) + x £p (< 7 , w) +x+(g,w) which is consistent with Eq. (6) of Ref. [18]. We em- 

s ---' phasize here that even though the density dependence 
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for massless and massive Dirac systems is different from 
that of the parabolic systems, they can be rewritten in 
the same form as Eq. (3), 

~ + o K). (io) 

£p /VfU 2 \UJ V 

Note the similarity between Eq. (7) and Eq. (10). This 
prompts the following mapping: band mass in parabolic 
systems, m p —> m a = £ F /vp, density dependent effective 
Dirac mass in massive as well as massless Dirac systems 
(to be distinguished from the band gap A which is also 
occasionally referred to as mass ). As a natural conse¬ 
quence, this correspondence manifests itself in all of the 
subsequent calculations. 

Physically the Dirac mass is a dynamical collective 
mass and is essential to explain inertial acceleration of 
the Dirac plasma under application of an external elec¬ 
tric field. In fact it has been insightfully defined as 
‘plasmon mass’ in the context of graphene 19 , and has 
also been recently measured in graphene 10 . Note that 
to d = £ f / v ^ is also the cyclotron effective mass (m c ) for 
Dirac systems 35 , which is typically defined as 27rm c = 
H 2 dS(£)/d£, where S(£ ) = 7r(e 2 — A 2 )/H 2 Vp denotes the 
area of a closed cyclotron orbit of massive Dirac electrons 
with energy £. For a system with a parabolic dispersion, 
the band mass, plasmon mass and the cyclotron mass are 
identical, m p = m c . 

Having obtained the long wavelength limit of the dy¬ 
namical polarization function, we now proceed to calcu¬ 
late the long wavelength limit of the plasmon dispersion 
in the next section. 


III. PLASMON DISPERSION 


Using Eq. (7) and Eq. (2) in Eq. (1), the well known 
long wavelength plasmon dispersion for systems with 
parabolic dispersion 2,5 , in one-, two- and three dimen¬ 
sions can be easily obtained to be 


^i p) = J 26 ni qV K o{qa) + 0(q 3 ) , 

V Mp 
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An important point to note here is that for the first term 
in Eq. (11), we can substitute to p —> to, i.e., replace the 
effective band mass by the classical mass of the particle 
and the quantum mechanical plasmon dispersion takes 
precisely the same form as that of classical density os¬ 
cillations in an electron liquid 1,2,5 . Physically this is a 
direct consequence of the fact that the long wavelength 
plasmons involve the motion of the entire plasma, and to 


lowest order it does not depend on the complex exchange 
and correlation effects that dress the motion of an indi¬ 
vidual electron. It should be emphasized that no such 
classical analog exists for Dirac systems, and the plas¬ 
mon dispersion in Dirac systems is intrinsically quantum 
mechanical in nature 18 . Note, however, that the higher 
order correction terms in Eq. (11) are fully quantum me¬ 
chanical and H explicitly appears in them. 

For systems with massive Dirac dispersion, using 
Eqs. (8) and (2) in Eq. (1), the q —» 0 limit of the plasmon 
dispersion is given by 
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+ 0{q 2 ), (12c) 


where we have defined g = g s g v . Note that 2d massive 
Dirac plasma was also studied in Ref. [24] , which reported 
an expression similar to Eq. (12b). For the limiting case 
of A —>• 0, Eq. (12) leads to, 


’i s) = 

\j htiij qV K 0 (qa) + 0(q 3 ) , 
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Equation (13) reproduces the results for gapless Dirac 
plasma reported in Ref. [18]. One important similarity 
between Eqs. (11), (12) and (13), is the same functional 
dependence of the plasmon frequency on the wave vector 
q. This is direct consequence of the physical requirement 
that the long wavelength plasmon dispersion must sat¬ 
isfy particle conservation (or continuity equation). One 
important difference between parabolic and Dirac sys¬ 
tems is that while the long wavelength limit of plasmon 
dispersion in parabolic systems is essentially ‘classical’ 
in nature, the plasmon dispersion in GDP and MDP is 
essentially quantum mechanical, as evidenced by the ex¬ 
plicit appearance of h in Eqs. (12) and (13). 

The long wavelength dependence of the dynamical po¬ 
larization function is the same for parabolic systems, 
GDP and MDP: n oc q 2 /oj 2 , however the proportion¬ 
ality constant has a different density dependence for 
various systems. For parabolic systems n( p ) oc 
for gapless Dirac systems n( g ) oc , and for mas¬ 

sive Dirac systems oc n.d/{n 2 J d + a^A 2 ) 1 / 2 , where 
= ( 9 A) 2 , g/ 4?r, (g/6n 2 ) 2/3 in 1, 2, and 3 dimen¬ 
sions respectively. As a consequence, the density depen¬ 
dence of the long wavelength plasmon frequency for MDP 
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is completely different compared to GDP and parabolic 
dispersion systems. As seen from Eq. (11), the plasmon 
frequency for parabolic dispersion system is proportional 
to y/rid in all dimensions. However, for GDP the plasmon 

dispersion follows aoc y/n d /n d ^ 2d behavior and for the 
one dimensional case, the plasmon mode is completely in¬ 
dependent of the density. For the MDP case, the density 
dependence completely changes due to the presence of 
gap and takes the form u>^ a \/^d/{ r k d d + ayA 2 ) 1 / 4 , 
where a d = (g/ 7r) 2 , g/ 47r, (g/6n 2 ) 2 / 3 , in 1, 2, and 3 di¬ 
mensions respectively. Note that in one dimension, the 
plasmon frequency in GDP is independent of the density 
whereas for MDP, the plasmon frequency of MDP has an 
explicit dependence on n\ as evident from Eq. (12a). 

Similar to the case of the long wavelength polarization 
function in Eq. (10), the plasmon frequencies in the q —> 0 
limit in GDP and MDP can also be rewritten in the form 
similar to Eq. (11) for parabolic systems, 


wj m,s) = J^-^qy/K 0 (qa) + 0(q 3 ) , 
y Krrid 

= j2ne 2 n l q 1 / 2 + 3/2 , 

y Km d 
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(14b) 
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We emphasize again that even though the density depen¬ 
dence of plasmon frequencies in long wavelength limit in 
the GDP and MDP is different from that of the parabolic 
systems, the frequencies can be rewritten in the same 
form using the density dependent effective Dirac mass (or 
equivalently the cyclotron mass). However, we note that 
this simplicity is deceptive, since for parabolic systems 
the band mass m p and density n d can also be treated 
as classical independent variables, but for Dirac systems 
m ( j is density dependent via the Fermi energy in a purely 
quantum mechanical way. 


IV. PLASMONS IN METAMATERIALS MADE 
OF MASSIVE DIRAC PLASMA, RIBBON AND 
LAYER, ARRAYS 


Coulomb interactions. Assuming no wave function over¬ 
lap between any consecutive layers or nano-ribbons, the 
collective modes of such superlattices, within RPA, are 
given by the zeros of the determinant of the general di¬ 
electric matrix of the superlattice, whose elements are 
given by 


ew = Su> - v w (<?, k)Ui' (q, w), (15) 

where n;(g,w) = n(q,w) is the bare density-density re¬ 
sponse function of each nano-ribbon or layer whose long 
wavelength limit is given in Eq. (8). In Eq. (15), vw(q , k) 
is the repulsive Coulomb interaction between the l and V 
nano-ribbon or layer in the periodic array which is given 

by 

vw = ^ ?re e~ qb ^ l ~ l I d = 2 , (16a) 

nq 

vw = — [SwK 0 (qa) + (1 - S w )K 0 (qb\l - l'\)] d = 1 , 

K 
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where b is the superlattice spacing. 

Assuming a periodic boundary conditions for the Id 
superlattice, the eigenvalues of the general dielectric ma¬ 
trix in Eq. (15) are given by 1 — v q Ii(q , uj)S d (q , k) where 
S d (q,k) = v~ 1 J2i' v W e ~ zk< ' l ~ l '^ b is the form factor for 
the Id superlattice formed from a d-dimensional plasma, 
and k is to be interpreted as a new wave vector arising 
form the periodicity of the infinite superlattice array and 
|fc| < 7 x/b. For definiteness we take the Id MDP wire 
array to be along the x axis ( q = q x ), and k = q y to 
be along the superlattice direction of the y axis. For the 
2d MDP layer superlattice, we consider it to lie in the 
x — y plane and the wave vector k = q z to be along the 
superlattice direction — the z axis. The dimensionless 
form factors can now be evaluated and are given by 


2 

-Si = 1 + ———r V I< 0 (nqb ) cos {nq y b) , (17a) 

M ^ 

a = Y' -q\n\b-iq z nb = sinh(qb) 

2 ^ cosh (qb) — cos (q z b) 


(17b) 


In this section we consider collective density excita¬ 
tions in metamaterials (periodic arrays) of massive Dirac 
plasma systems. In particular we consider a stacking 
of identical Id massive Dirac plasma nano-ribbons (or 
quantum wires) placed parallel to each other in a plane, 
and a periodic array of parallel 2d massive Dirac plasma 
sheets. Similar systems made of parabolic dispersion 
and gapless Dirac plasma have been theoretically investi¬ 
gated earlier 18 ’ 28 ^ 30 ’ 33 and experimentally demonstrated 
for graphene 8,9,11 . To describe the collective modes of 
such superstructures, we need to include the inter-ribbon 
or inter-layer Coulomb interactions, which leads to a cou¬ 
pling of all the layers due to the long range nature of 


Note that for q z = 0, as q —> 0, S 2 —> q- The plas¬ 
mon bands for the Id superlattice, composed of 

Schrodinger electrons, MDP and GDP, are now explic¬ 
itly given by the zeros of the eigenvalues of the general 
dielectric matrix in Eq. (15): 

1 - v q Il(q,uj)S d (q, k) = 0 . (18) 

However in the long wavelength dynamical limit {q — > 0), 
n oc q 2 /uj 2 for parabolic, massive Dirac and gapless Dirac 
systems in all dimensions, and consequently Eq. (18) sim¬ 
plifies to give, 


(19) 
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FIG. 1: (a) The long wavelength plasmon dispersion for 

Id and 2d MDP against the backdrop of the Id particle-hole 
continuum (the grey shaded region). The green shaded region 
marks the plasmon band formed in the Id superlattice of the 
Id MDP [the lower boundary is for q y = 7 r/6, and the upper 
boundary is for q y = 0 in Eq. (19)]. (b) The long wavelength 
plasmon dispersion for 2d and 3d MDP, with the green shaded 
region marking the plasmon band formed in the Id superlat¬ 
tice of the 2d MDP [the lower boundary is for q z = vr/fo, and 
the upper boundary is for q z = 0]. The gray shaded region 
in panel (b) depicts the particle-hole continuum of 2d and 
3d massive plasma. Note that as in the parabolic case, the 
electron-hole continuum in Id differs from the 2d and 3d case 
because there are no excitations at finite q and low u> even 
for massive and massless Dirac fermions. In both panels we 
have defined hvpko = ep and used the following parameters: 
A/sf = 0.4, ge 2 / (2Hkvf) = 2, ako = 0.25, and bko = 2. 


where the form factor Sd is explicitly given in Eq. (17). 
We emphasize that Eq. (19) is very general and it de¬ 
scribes the long wavelength plasmon dispersion for Id 
superlattice made of parabolic, massive Dirac or gapless 
Dirac systems (for both ribbons and layers). The plas¬ 
mon bands for a superlattice of Id MDP nano-ribbons, 
and 2d MDP layers is displayed in panels (a) and (b) of 
Fig. 1, respectively against the backdrop of the particle- 
hole continuum. 

For ribbons and layers of parabolic systems, the super¬ 
lattice plasmon dispersion at the upper band edge (k = 0 ) 


is given by 


q y = °) 

^2s ( 9 ; Qz = 0 ) 


/ 2'KTl2e 2 q 
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( 20 a) 

( 20 b) 


with h 2 = qy and 713 = r y. Note that in Eq. (20), the d- 
dimensional superlattice plasmon at the band edge (fc = 
0 ) has exactly the same form as the corresponding (d+ 1 )- 
dimensional bulk plasmon [see Eqs. (lla)-(llb)] with the 
effective densities being n 3 = n 2 /b and n 2 = ni/b. This is 
consistent physically since the d-dinrensional superlattice 
loses its discrete periodicity at the band edge (k = 0 ) and 
effectively becomes a (d+l)-dimensional system. 

In the case of superlattice structures made of MDP, 
the q —> 0 plasmon dispersion at the band edge (k = 0 ) 
is 


'l?\q\Qy = 0) = 1 

l2ge 2 VFq (e| — A 2 ) 1 / 4 
v hftb p 1 / 2 

b F 

( 21 a) 

’ 2 ^ ( 9 ; 9* = 0) = 1 

1 gc 2 Is 2 - A 2 

( 21 b) 
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These in the limit A —> 0, give the corresponding expres¬ 
sions for the GDP: 

/ 2 \ 1/2 

^W = 0 )=(^) . ( 22 b) 


The physically appealing correspondence between the d- 
dimensional superlattice at the band edge and (d+ 1 )- 
dimensional bulk system, does not happen for MDP 
as well as for GDP. From Eqs. (21a)-(21b), it is clear 
that at the band edge, uj[™ ,9 \q;q y = 0 ) 7 ^ with 

the intuitive substitution h 2 = n\/b and u>^’ 9 \q-,q z = 
0) 7 ^ u^ n ' 9 ' > with fi 3 = n 2 /b. This is a direct conse¬ 
quence of different density dependence of the polarization 
function [see Eqs. ( 8 ) and (9)] in massive and massless 
Dirac systems as compared to systems with a parabolic 
dispersion relation [see Eq. (7)]. However, the 2d su¬ 
perlattice plasmon dispersion would map to the corre¬ 
sponding plasmon dispersion for 3d massive Dirac plasma 
u> 2 ^\q;q z = 0 , 712 ) —> £*4”^ (q,n 3 ), if rather than the in¬ 
tuitive definition h 3 = n 2 /b , we have the following corre¬ 
spondence of densities in 3d and 2d massive systems, 

f 2 = (67r 2 n 3 g- 1 ) 2/3 + A 2 

\ n 2 ) + A 2 

In the A —> 0 limit, Eq. [23] leads to the corresponding 
relation GDP, i.e., h 3 = (9Trg/16) 1 ^ 4 (n 2 /b 2 ) 3 ^ 4: , which 
was first derived in Ref. [18]. For the Id superlattice, the 
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gapless Dirac plasma frequency w[ 9 J does not depend on 
the carrier density at all, and hence massive Dirac plasma 
frequency differs here in this aspect. For Id massive Dirac 
plasma, the superlattice plasmon frequency at the band 
edge would agree with the 2d massive plasma frequency, 
i.e., iu[™\q;q y = 0,ni) —> w?™'* (q,ni), only if we have 
the following relation between the densities fi 2 and n\. 

f h 2 b\ 2 = 4t m 2 g + g 2 A 2 

V n i ) TT 2 nf+g 2 A 2 


In the A —>■ 0 limit, Eq. [24] reproduces the corresponding 
relation GDP, i.e., ro 2 = Ag/nb 2 . 

Due to the presence of gap in massive Dirac systems, 
the density dependence as well as the band edge plas¬ 
mon at k = 0 differs from the usual parabolic as well as 
the massless Dirac systems. This is because of the fact 
that the density dependence of the polarizability of mas¬ 
sive Dirac systems is completely different as compared 
to massless Dirac and parabolic dispersion systems, as 
shown in Eqs. (8). 

Finally, we note that as in the case of bulk plasmons, in 
Eqs. (14a)-(14c), the superlattice plasmon frequencies for 
Dirac systems can be expressed in the same form as for 
systems with parabolic dispersion relation. Expressing 
the numerator in Eqs. (21a)-(25a) in terms of density, 
and the denominator in terms of the cyclotron mass of 
massive Dirac particles = ep/up, we have 

(m, g ), [2irn 2 e 2 q \ 1/2 

-L w = °) = (^ mr j - ( 25a ) 

w 2 7’ 9 \^Qz = 0) = ( 4lrri3e ) . (25b) 

V ™ d ) 


V. SUMMARY AND CONCLUSION 

In this article we have obtained the long wavelength 
plasmon frequency for massive Dirac particles in various 
dimensions and their Id superlattice, and compared it 
with the corresponding results for plasmons in parabolic 
systems and gapless Dirac systems. As expected, fac¬ 
tors of 1/y/h explicitly appear even in the leading order 
term in the long wavelength plasmon dispersion of MDP 
and GDP, highlighting their intrinsically nonclassical and 
quantum nature. 

To summarize we find that the long wavelength limit 
of the dynamical density response function, while hav¬ 
ing the same dependence for q and w, i.e., II oc q 2 /uj 2 


for parabolic as well as Dirac systems, has a different 
dependence on density. This different density depen¬ 
dence also gets manifested in the long wavelength plas¬ 
mon dispersion and we find that for massive Dirac sys¬ 
tems ui^ oc \fndl{n'J L d + o^A 2 ) 1 / 4 while for gapless 
Dirac systems oc y/nd/‘ri 1 J 2d , and for parabolic sys¬ 
tems oc y/rid, in d-dimensional systems. Addition¬ 
ally we note that a beautiful similarity emerges between 
all three systems if we use the density dependent effec¬ 
tive Dirac mass (or cyclotron mass) for GDP and MDP, 
ro<j = ep/r’l, to express the long wavelength plasmon 
dispersion for all systems in all dimensions: a/ p ’ s,n d oc 
y/rid/m p/d- This density dependence of the plasmon fre¬ 
quency may be used to distinguish between various types 
of systems (parabolic, GDP and MDP) and their effective 
dimensionality. Alternatively, the long wavelength plas¬ 
mon dispersion may be used to determine the dynamical 
collective mass of various Dirac systems 10 . 

We have also calculated the plasmon dispersion in ald- 
superlattice made of nano-ribbons and layers of MDP, 
and find that while the q dependence is similar to that 
of parabolic superlattice, the density dependence is com¬ 
pletely different. Note that our results for the 3d plas¬ 
mons of MDP, 2d layers and multilayers, as well as Id 
ribbons and multi ribbon arrays can be tested using elec¬ 
tron scattering, light scattering or infrared spectroscopy 
for ribbons and layers made of silicene, transition metal 
dichalcogenides and other materials which have a disper¬ 
sion similar to that of massive or gapped Dirac spectrum 
at low energies. 

Finally we note that we have not included in our calcu¬ 
lation the effect of dielectric mismatch, which may lead 
to the dielectric constant appearing in the expression of 
the Coulomb interaction k, to be wave vector dependent. 
This can be easily included in our calculations according 
to Ref. [36] . However we believe that the long wavelength 
plasmon dispersion, as discussed in this article, would not 
be impacted by this, even though the dielectric mismatch 
may change the plasmon dispersion for finite wave vec¬ 
tors. 
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